/***********************************************************
Compile: gcc -static -I. -L. -O3 grn.c -o grn -liode_intel64 -lm
Usage  : ./grn N input (./grn 10 input)

Readme : Solve ODE for Gene Regulatory Network
***********************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "intel_ode.h"

double Gamma=0.1;
double u=0.0;
double h=2.0;
double K=3;
double t_step=1.e0;

int *network;

void grn(int*n,double*t,double*y,double*f)
{
	int i, j, N;
	double A, R;

	N = *n;
	for (i=0;i<N;i++)
	{
		A = R = 0.0;
		for (j=0;j<N;j++)
		{
			if (i==j)
				continue;
			if (network[i*N+j] == 1)
				A += y[j];
			else if (network[i*N+j] == -1)
				R += y[j];
		}
		A = pow(A,h)/(pow(A,h)+pow(K,h));
		R = pow(K,h)/(pow(R,h)+pow(K,h));
		f[i] = u-Gamma*y[i]+(1.0-u)*A*R;
	}
}

void solve_ODE (int N, int init_state[N])
{
	int n, ierr, i, j;
	int kd[N], ipar[128];
	double t, t_end, h, hm, ep, tr;
	double y[N], dpar[(7+2*N)*N];

	hm=1.e-12; /* minimal step size for the methods */
	ep=1.e-6;  /* relative tolerance. The code cannot guarantee
		      the requested accuracy for ep<1.d-9 */
	tr=1.e-3;  /* absolute tolerance */

	for (i=0;i<128;i++) ipar[i]=0;
	t=0.e0;
	h=1.e-7;

	n=N;
	t_end=t_step;
	for (i=0;i<N;i++) y[i]=init_state[i]/Gamma;

	for (i=0;i<100;i++)
	{
		dodesol_rkm9mkn(ipar,&n,&t,&t_end,y,grn,&h,&hm,&ep,&tr,dpar,kd,&ierr);

		if(ierr!=0)
		{
			printf("dodesol_rkm9mkn routine exited with error code %4d\n",ierr);
			exit (-1);
		}

		printf ("t=%5.1lf", t_end);
		for (j=0;j<N;j++)
		{
			printf (" %5.3lf", y[j]);
		}
		printf ("\n");

		for (j=0;j<128;j++) ipar[j]=0;
		h=1.e-7;
		n=N;
		t=t_end;
		t_end+=t_step;
	}
}

int main (int argc, char *argv[])
{
	FILE *fp;
	char buf[256];
	int i, j, N;

	N = atoi (argv[1]);

	int init_state[N];
	network = calloc (N * N, sizeof (int));

	fp = fopen (argv[2], "r");
	while (fgets (buf, 256, fp))
	{
		for (i = 0; i < N; i++)
		{
			fscanf (fp, "%d", &init_state[i]);
		}
		fgets (buf, 256, fp);

		for (i = 0; i < N; i++)
		{
			fgets (buf, 256, fp);
			for (j = 0; j < N; j++)
			{
				if (buf[2 * j] == 'n')
					network[i * N + j] = 0;
				else if (buf[2 * j] == 'g')
					network[i * N + j] = 1;
				else if (buf[2 * j] == 'r')
					network[i * N + j] = -1;
				else
				{
					printf ("Error in reading network topology!\n");
					exit (-1);
				}
			}
		}

		solve_ODE (N, init_state);
	}
	fclose (fp);
}
